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FAR  FIELD  NUMERICAL  BOUNDARY  CONDITIONS  FOR  INTERNAL 
AND  CASCADE  FLOW  COMPUTATIONS 

Ch.  HIRSCH, 

Department  of  Fluid  Mechanics, 

Vrije  Universiteit  Brussel,  Belgium 


ABSTRACT 

The  present  report  extends  the  approach  developed  by  A.  Verhoft  tor  the  treatment  of  the 
far  field  boundary  conditions,  Verhoff  and  O'Neil  (1984),  to  more  general  formulations  of  the 
Euler  equations  and  to  cascade  geometries. 

Linearized  solutions  of  the  Euler  equations  are  developed  for  the  perturbations  from  the 
uniform  free  stream,  for  ducts  and  cascades.  These  solutions  are  based  on  the  conditions  that 
the  waves  associated  with  incoming  characteristics  should  decay  to  zero  in  the  far  field,  while  the 
variables  associated  to  the  outgoing  characteristics  are  derived  from  the  numerical  internal 
solution.  The  exact  linearized  solutions  are  based  on  a  Fourier  expansion  in  the  direction  along 
the  inlet  or  exit  boundaries. 

Results,  obtained  from  an  Euler  code  are  shown  for  ducts  and  cascades,  comparing 
calcualtions  for  exit  boundaries  at  increasingly  closer  distance  to  the  central  flow  region.  The 
method  is  also  valid  for  nonisentropic  flows. 

The  results  show  that  the  corrections  to  the  uniform  boundary  conditions  derived  from  the 
analysis  allow  a  considerable  reduction  of  the  computational  domain,  with  the  corresponding 
savings  in  computational  times. 


FAR  FIELD  NUMERICAL  BOUNDARY  CONDITIONS  FOR  INTERNAL 
AND  CASCADE  FLOW  COMPUTATIONS 


Ch.  IIIRSCII, 

Department  of  Fluid  Mechanics, 
Vrije  Universiteit  Brussel,  Belgium 


INTRODUCTION 

The  present  report  extends  the  approach  developed  by  A.  Verhoff  for  the  treatment  of 
the  far  field  boundary  conditions,  Verhoff  and  O'Neil,  (1984),  Verhoff  (1988),  to  more 
general  formulations  of  the  Euler  equations  and  to  cascade  geometries. 

The  imposition  of  uniform  boundary  conditions  in  the  subsonic  far  field  of  the 
computational  domain,  such  as  constant  pressure  condition,  requires  the  computational 
boundary  to  be  located  at  a  large  distance  from  the  central  flow  region.  Although  the  flow 
variations  in  the  vicinity  of  the  external  boundaries  are  generally  smooth,  a  large  number 
of  mesh  point  are  necessary  to  cover  adequately  these  extended  regions. 

If  more  information  would  be  available  on  the  behavior  of  the  flow  in  the  far  field, 
allowing  for  non-uniform  flow  distributions  along  boundaries,  it  would  be  possible  to 
accept  external  boundaries  at  distances  much  closer  to  the  main  flow  region,  with  the 
subsequent  reduction  in  number  of  mesh  points. 

In  this  report  linearized  solutions  of  the  Euler  equations  are  developed  for  the 
perturbations  from  the  uniform  free  stream,  for  ducts  and  cascades.  These  solutions  are 
based  on  the  conditions  that  the  waves  associated  with  incoming  characteristics  should 
decay  to  zero  in  the  far  field,  while  the  variables  associated  to  the  outgoing  characteristics 
are  derived  from  the  numerical  internal  solution.  The  exact  linearized  solutions  are  based 
on  a  Fourier  expansion  in  the  direction  along  the  inlet  or  exit  boundaries. 

Typical  outcome  of  the  analysis  leads,  for  a  cascade  exit  station  with  subsonic  axial 
velocities,  to  the  replacement  of  the  constant  pressure  condition  by  an  improved  boundary 
condition,  relating  the  pressure  perturbations  to  the  Fourier  expansion  of  the  flow  angle 
along  the  exit  boundary  calculated  from  the  numerical  internal  solution. 

Results  are  also  presented  for  ducts  and  cascades,  comparing  the  results  for  exit 
boundaries  at  increasingly  closer  distance  to  the  cascade  blades. 

In  the  first  part,  sections  1  and  2,  we  present  a  reformulation  of  A.  Verhoff  s  (1985) 
isentropic  analysis  for  ducts,  based  on  the  perturbation  analysis  of  the  general  form  of  the 
characteristic  equations,  as  applied  by  llirsch  et  al.(  1987).  The  first  order  equations  are 
identical  to  those  derived  by  Verhoff  for  his  generalized  Riemann  invariants,  although  the 
acoustic  wave  perturbations  are  represented  by  a  different  combination  of  variables.  This 
covers  therefore  isentropic  as  well  as  non-isentropic  flows.  As  a  by-product,  an 
interesting  relation  is  derived,  relating  the  far  field  pressure  and  velocity  fluctuations  for 
isentropic  flows. 

The  third  section  presents  the  extension  of  the  approach  to  cascades.  Due  to  the 
periodicity  boundary  conditions  in  the  inlet  and  outlet  domains,  resp.  upstream  and 
downstream  of  the  cascade  blades,  a  complex  Fourier  representation  has  to  be  defined  for 
the  variations  in  the  direction  of  the  cascade  front.  Results,  obtained  from  an  Euler  code 
are  shown  for  ducts  and  cascades,  confirming  the  correctness  of  the  analysis  and  the 
potential  CPU  benefits  of  closer  boundaries. 


1. PERTURBATION  ANALYSIS  FOR  GENERALIZED 
CHARACTERISTIC  EQUATIONS. 

The  small  perturbation  formulation  of  the  Euler  equations  can  be  derived  in  many 
ways,  depending  on  the  particular  variables  considered  for  the  analysis.  The  most  current 


formulation  puts  forward  the  pressure,  since  it  is  the  basis  for  all  acoustic  applications, 
and  takes  the  Euler  equations  in  primitive  variables,  p,p  and  velocity,  as  starting  point. 

In  computational  fluid  dynamics,  appropriate  information  on  outgoing  pressure 
radiation  waves  allows  to  investigate  far  field  boundary  conditions,  which  should  prevent 
spurious  waves  to  be  reflected  from  the  boundaries  towards  the  computational  domain. 
This  has  been  investigated  by  Bayliss  and  Turkel  (1982),  where  references  can  also  be 
found  to  earlier  important  work  in  that  direction  and  more  recently  by  P.Roe  (1986).  In 
these  methods,  differential  equations  for  the  pressure  are  derived  for  the  far  field  pressure 
and  added,  after  appropriate  discretization,  to  the  other  internal  equations. 

Other  approaches  make  use  of  the  far  field  formulation  of  the  small  perturbation 
potential  equation  and  introduce  analytical  asymptotic  expansions  as  boundary 
corrections,  Thomas  and  Salas  (1986). 

In  the  present  approach,  numerical  information  is  used  directly  in  interaction  with  a 
series  expansion  of  the  far  field,  applying  a  perturbation  expansion  of  the  characteristic 
equations.  Hence,  we  take  as  starting  point,  the  compatibility  equations.  Note  that  these 
different  point  of  views  are  clearly  interconnected:  Roe  (1986)  derives  a  particular  form  of 
the  acoustic  compatibility  equation  corresponding  to  a  selected  wave  propagation  direction 
which  corresponds  to  the  selected  pressure  equation.  In  the  present  work,  far-field 
pressure  relations  will  be  derived  from  the  characteristic  analysis. 


1.1.  Euler  equations  in  characteristic  form 

We  consider  the  general  form  of  the  Euler  equations  in  characteristic  form, 
written  for  an  arbitrary  direction  of  wave  propagation  ic>  considered  as  a  unit  vector, 
following  Hirsch  et  al.  (1987),  see  also  Hirsch  (1988),  (1989). 

The  dependent  variables  are  the  following  characteristic  quantities,  defined  in 
differential  form, 


aw,  =  ap  ---dp 

(1.1a) 

0W2  =  l.av  =  Ky0U  -  Ky0V 

(1.1b) 

0W  3  —  ic.av  +  “dp=  K/)U  +  K  y0V  + 

(1.1c) 

dw4  =  -  K.0V  +  ^,0p=  -  (tCjflll  +  K^V)  +  p~0p 

(lid) 

where  the  vector  I  is  normal  to  the  wave  number  vector  ic  and 
tcx)  r  Note  that  1  is  also  a  dnit  vector;  see  figure  1.1. 

has  the  components  (Ky 

The  first  characteristic  variable  is  proportional  to  the  entropy  and  the  associated 
equation  will  describe  the  convection  (propagation)  of  entropy  (waves).  The  second 
component  is  the  amplitude  of  a  vorticity  or  shear  wave,  which  has  no  equivalent  in  one 
dimensional  flows.  The  third  and  fourth  components  are  the  amplitudes  of  the  acoustic 
waves. 

The  compatibility  equations  for  the  Euler  system  of  conservation  laws,  become 


0,w  |  +  (v.^Ow!  =  0 


(1.2a) 


T 


(1.2b) 


3,w 2+  (v.V)w2  +  J  (I.VXW3+W4)  =0 
dtw3+  (v+cjc).^w3+  c  (l.^)w2  =  0  (l-2c) 


atw4  +  (v-ctc).^w4+  c  (l.?)w2  =  0  (1 ,2d) 

The  first  two  terms  of  all  the  equations  (1.2)  are  purely  convective  and  represent  the 
propagation  of  the  associated  wave  in  the  characteristic  directions.  The  third  terms 
represent  the  coupling  or  the  interaction  between  the  different  waves  and  their  presence 
results  from  the  fact  that  the  jacobian  matrices  afe  not  simultaneously  diagonalizable  for 

wave-like  solutions  of  the  formU=Uo  e  x  elt01  whereby  all  the  components  of  U,  that 

is  all  the  flow  variables,  propagate  along  the  same  spatial  direction  k. 

If  this  direction  is  taken  aligned  with  the  velocity,  than  the  above  equations  can  be 
transformed  to  the  particular  streamwise  characteristic  form  used  by  Verhoff,  Verhoff  and 
0'NeiI(1984).  Defining  the  dimensionless  entropy  as 

S  =  -s/(yr)  (1.3a) 

where  s  is  the  physical  entropy  and  r  the  gas  constant,  we  have 


S-S  = 


_  -1 


Y(Y-1> 


log 


P/P„ 


M 


with 


dS 


P(Y-1) 


dw  1 


(1.3b) 


(1.3c) 


The  constant  is  taken  as  =2/(y-l). 

Hence,  the  first  characteristic  equation  reduces  to  the  convection  of  entropy, 

as  as  / 1 

^  +  =  (,'4) 

where  ds  is  the  elementary  arc  length  along  the  flow  pathline  and  q  the  magnitude  of  the 
velocity. 

The  second  characteristic  variable 

dw2=ldv  (1.5) 

reduces  to  the  following  form  for  a  direction  K  -  I  v  ,  the  unit  vector  in  the  direction  of  the 

velocity  and  I  =  1  n  the  unit  vector  in  the  direction  nomial  to  the  velocity,  figure  1.2. 

With 


dv  =  dqT v+qd8T  n 


(1.6) 


9w2  =  1y.3v 


(1.7) 


the  second  equation  (1.2b)  becomes 

f  +  =  °  d.8) 

where  P  is  the  logarithm  of  the  pressure  and  9/9n  is  the  derivative  in  the  direction  normal 
to  the  velocity. 

The  third  characteristic  variable  is  written  in  the  following  form 

dw3  =  j^  + Vdv  =  “  +  aq  (1.9) 

and  equation  (1.2c)  becomes 

9w3  ,  9w3  09 

_aT  +  (q  +c)Ts^  +qcaTi =  0  ( 1  • 1 0) 


Introducing  equation  (1.9)  and  expressing  the  pressure  variations  dp  in  function  of  the 

variations  of  entropy  and  speed  of  sound,  via  the  definition  c2=rp/p  for  a  perfect  gas  and 
equation  (1.3),  leads  to 


~  =  cdS  +  ---dc  =  d(cS)  -  (S  -  —  )dc 
pc  Y-l  '  v  Y-l 


(1.11) 


The  formulation  of  the  Euler  equations  used  by  Verhoff  is  based  on  generalized' 
Riemann  variables  R  and  Q,  defined  as  R=q-cS  and  Q=q+cS  and  equation  (1.10)  can  be 
transformed  to 


9(q+cS) 


9(q+cS) 


2  Jc 


,9c,  99 


+  (q+c)^ir-  -  (s  -  ♦  (q+c>i)  - 


at  w  9s  Y-rl3t 

Note  also  the  interesting  telation  for  the  speed  of  sound 

9c  9c  Y-1 

9t-  +  q9s  =  -"2C(Vv) 


(1.12) 


(1.13) 


and  with  the  expression  of  the  divergence  of  the  velocity  in  local  coordinates, 


3  -  9q  90 
v.v  =  ^.-  +  q=r 
9s  M9n 


(114) 


equation  (1.12)  reduces  to  the  form 


9(q+cS) 


+  (q+c)- 


9(q+cS) 


/c  2  ,  ,9c  Y-l  9q  Y-l  o90 
=  S  y_  |  C  5s  2  .«> '  2  <ioS.n 


(1.151 


used  by  Verhoff  and  0'NeiI(1984). 

The  fourth  characteristic  equation  leads  to  the  similar  equation  for  the  variable  R=q-cS 


0(q'cS>  i  (a  cl8{q'cS)  t*'  2  )  c  (9c  ,  Y_I  dqi  ,  y-lac*dQ 

at  +(q'c)  as  -{S_Y_1)c(as+  2  as)+  2  qcSan 


(1.16) 


Equations  (1.4),  (1.8),  (1.15)  and  (1.16)  are  identical  to  the  'Quasi  ID’  formulation  of 
Verhoff  and 
O'Neil  (1984). 


1.2.  Perturbation  expansion  in  the  far  field 

If  the  flow  at  infinity  is  uniform,  at  any  finite  distance  a  perturbation  will  exist,  which 
is  supposed  to  be  small  at  some  distance  from  the  central  region  of  the  flow  domain. 
Denoting  by  a  subscript  the  uniform  quantities  and  writing  all  variables  U  as 

U=UOQ+eU’,  where  e  is  a  small  perturbation  parameter,  the  first  order  linearized  Euler 
equations  (1.2)  take  the  form 


djuv,  +  (V00V)w1  =  0 

(1.17a) 

atw2  +  (\L.v)w2  +  ~  (I.v)p'  =  0 

(1.17b) 

a,w3  +  (V00+  cmk).Vw3  +  cJ.V)w2  =  0 

(1.17c) 

a,w4+  (V„-c00ic).Vw4  +  c00(i.V)w2=  0 

(1 •  17d) 

The  perturbation  equations  are  written  here  in  the  general  form  for  an  arbitrary  wave 
propagation  direction  and  the  first  order  corrections  w’  of  the  characteristic  variables  are 
defined  as  follows,  in  function  of  the  perturbations  S',  p'  and  v'  of  resp.  entropy, 
pressure  and  velocity.  The  first  characteristic  variable  w'j  is  proportional  to  the  entropy, 
hence  we  have,  following  equation  (1.3b), 

w'j  =  p(  y*  1  )S'  (1.18) 

while  the  second  characteristic  variable  w'2.  in  accordance  with  equation  (1.5),  is  defined 
as 


w2  =  l.v'  =  v'cosa  -  u'sina  (1.19a) 

referring  to  figure  1.3,  a  being  the  angle  bctwen  the  propagation  direction  k  and  the 
reference  x-direction  and  u',  v'  the  cartesian  components  of  the  perturbation  velocity.  In 
function  of  the  flow  angle  deviation  0'  to  the  free  stream  angle  8M,  and  the  projection  q 


of  the  velocity  perturbation  on  the  direction  of  the  far  field  velocity,  the  above  equation 
can  also  be  written  as 


w2=  ^e'cosfa-ej-q's^a-ej  (1.1%) 

The  perturbations  of  the  third  and  fourth  characteristic  variables  are  defined  by 

w3  =  -£—-  +  qooe’  sin(a-0cJ  +  q'cos(a-eoo)  (1.20) 

P  oo 

*  oo 


w4  =  -^r--qoo9'  sin(a-eoo)-q'cos(a-eoo)  (1.21) 

P«~> 

Transforming  the  time  variable  to  a  variable  with  space  dimensions,  via 

T  =  Coot  (1.22) 

and  introducing  the  free  stream  Mach  number 

M 

OO  =  q°o/c  oo  (1.23) 

leads  to  the  following  form  of  the  first  order  perturbation  equations,  where  x  and  y  are 
reference  cartesian  coordinates. 

j^  +  MJcosO^  +  sind^lS^O  (1.24a) 

+  M  JcosO„,^-  +  sin0„^)w2-  -  (sina  --  -cos«  ^)(w3+w4)  =  0  ( 1  24b) 

dw3  dw3  i)w3  r)w2  dw? 

+  (M00cos000+  cosa)—  4  (M^sinO^  +  sina)  —  -  sina  -  +  cost* -j~  =  0  (1  24c) 


0W4  <lw4  dw.  Dw?  r)w2 

+  (M„,cos0„-cosa)  4  (M^sinO^-sinn)— --  sina  —  +  cos«  ^ --  =  0  (1.24d) 

Note  that  the  average  of  the  acoustic  wave  perturbations  is  proportional  to  the  pressure 
disturbance  p' 


(w3+  w4 


)/2  = 


(1.25) 


and  that  the  second  equation  (1.24b)  relates  the  flow  angle  perturbations  9'  to  the  far  field 
pressure  disturbance. 


A  simplified  form  is  obtained  when  the  wave  propagation  direction  ic  is  aligned  with 
the  free  stream  velocity  direction,  that  is  for  a=0oo.  In  addition,  the  x-direction  can  also 
be  aligned  with  the  free  stream  velocity,  that  is  a=9oo=0  and  the  above  equations  simplify 
considerably. 

The  characteristic  variables  reduce  to 


W,  =  p(y-1)S' 

(1.26a) 

w2=  q„e* 

(1.26b) 

w3  =  -£ —  +  u' 

(1.26c) 

3  p~c°° 

'  p' 

w  .  =  — - - u 

(1.26d) 

4  P„C°o 

and  the  system  (1.24)  becomes 

3S'  .  3S'  .  . 

+  (1-27a) 


8w2 


t  I  I 

3w2  j  3(w3+w4) 

ir  +  2"~^r“ 


=  0 


(1.27b) 


9w3 

~5T 


+ 


(Mm+1) 


3w3  3w2 

+  W 


=  0 


(1.27c) 


8w4 

3t 


+ 


3w,  3w. 


(m00-i)-3-2  + 


3x  3y 


(!  27d) 


This  formulation  of  the  first  order  perturbation  equations  is  identical  to  the  isentropic 
formulation  applied  by  Verhoff  (1985). 


1.2.2.  Relation  between  pressure  and  velocity  disturbances 

An  interesting  relation  connecting  pressure,  velocity  and  entropy  disturbances  can  be 
obtained,  to  first  order,  from  energy  conservation  or  constancy  of  stagnation  enthalpy. 

Writing  h+q2/2  =  llM  and  expanding  to  first  order  leads  to  the  following  relation 
between  enthalpy  and  velocity  perturbations 

h'  +  q'qso  =  0  (1.28) 
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where  q'  is  the  perturbation  of  the  magnitude  of  the  far  field  velocity.  This  quantity  is  not 
equal  to  the  magnitude  of  the  velocity  disturbance  v '  and  we  have  q'~u'. 

Expanding  the  definition  of  enthalpy  in  function  of  entropy  and  pressure 


Y-l 
_h 
h0 

leads  to 

s; 

Y  P~+CP 


(1.29) 


(1.30) 


Combining  with  equation  (1.28),  we  obtain  with  the  perfect  gas  relation  h00=  c^,/tY-l) 
and  the  dimensionless  entropy  disturbance  S', 


P  Cc 


+  q'Moo-cooS’  =  0 


(1.31) 


The  interest  of  this  relation  is  that  it  expresses  the  far  field  influences  of  entropy 
fluctuations  on  the  relations  between  the  velocity  and  pressure  fluctuations. 

The  acoustic  wave  amplitudes  dwj  and  dw4  can  also  be  expressed  in  function  of 
entropy,  speed  of  sound  and  velocity,  after  elimination  of  the  pressure  through  the 
entropy  relation  (1.11),  leading  to 

dw3  =  cdS  +  ^-dc  +  dq  (1.32) 

if  the  wave  propagation  direction  ic  is  aligned  with  the  local  velocity.  Similarly, 

dw4  =  cdS  +  ~ydc  -  dq  (1.33) 


The  perturbations  of  these  quantities  can  readily  be  defined  and  also  can  be  combined 
with  equation  (1.31),  to  give 


w„  =  — 
3  P. 


Or' 

+  a' =  c  S'  +  ---  +  a' 

■  oo  Y  1  ' 


=  cooS'  +  (1-MJq’ 


and 


w,  =  — - q'  =  cMS'  +  ----- -  q‘ 

4  P„cOT  M  ”  Y-l  M 


(1.34) 


=  C«,S  -  (1 +M  Jq’ 


(1.35) 


These  equations  imply  the  following  relation  between  velocity  and  speed  of  sound 
perturbations, 


2c: 

7-1 


(1.36) 


which  is  obtained  from  equation  (1.28)  by  an  expansion  of  the  perfect  gas  relation 

h=c2/(y-1)  (1.37) 

writing 

C  =  Coo  +  c’  (1.38) 


2.  FAR  FIELD  BOUNDARY  CONDITIONS  FOR  INTERNAL  DUCT 
FLOWS 

The  determination  of  far  field  boundary  conditions,  based  on  the  expansion  of  the 
characteristics  small  perturbations,  is  described  in  this  section  for  duct  flows,  with 
subsonic  in-  and  outflow  conditions. 

The  derivation  follows  the  approach  of  Verhoff  (1985)  and  is  repeated  here  with  the 
w'  variables,  leading  to  some  interesting  relations  for  the  pressure  and  velocity 
corrections. 

The  method  determines  analytical,  steady  state  solutions  for  the  perturbation  wave 
amplitudes  with  exponential  decay  in  the  streamwise  direction  (taken  as  the  x-direction) 
and  a  Fourier  expansion  in  the  direction  of  the  boundary  surfaces  (taken  as  the  y- 
direction).  The  unknown  coefficients  of  the  solutions  are  determined,  on  one  hand,  by 
expressing  that  the  incoming  characteristic  disturbances  at  the  boundaries  decay 
exponentially  in  the  far  field  and  on  the  other  hand,  by  matching  the  analytical  solution  at 
the  boundary  of  the  computational  domain  with  the  numerical  solution  for  the  outgoing 
characteristics. 

The  analysis  is  performed  for  stationary  flows,  for  the  simplified  formulation  (1.27) 
and  a  duct  geometry  represented  in  figure  2.1. 

The  duct  extends  to  infinity  with  a  computational  region  limited  by  inlet  and  exit 
boundaries  AB  and  CD,  the  width  of  the  channel  being  denoted  by  b.  We  look  for 
perturbation  solutions  of  the  system  (1.27)  in  the  region  between  the  inlet  (or 
exit)  station  of  the  computational  domain  and  the  boundary  at  infinity. 

Since  entropy  is  purely  convected  and  decoupled  from  the  other  equations,  we  can 
solve  separately  for  the  entropy  perturbation  and  remove  the  corresponding  equation  from 
the  system  (1 .27). 

Disturbance  solutions  by  separation  of  variables  are  sought  for  the  remaining  variables 
in  the  external  region,  with  a  Fourier  scries  in  y,  of  the  form 


w2=  X  hjx)  sin^  (2.1a) 

m«1 
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w3  =  X  U*)  cos^  (2- lb) 

m-1 


w4=  Egm(x)cosn?  (2-,c) 

m-1 

The  choice  of  the  Fourier  terms  in  equation  (2.1a)  results  from  the  flow  tangency 
boundary  condition  at  the  solid  walls,  requiring  that  0’=()  for  y  =  ±b/2.  Introducing  these 
solutions  in  the  stationary  form  of  equations  (1.27)  leads  to  the  following  system,  for  an 
arbitrary  Fourier  mode  m,  writing  M  instead  of  Moo,  and  removing  the  subscript  m  on  the 
amplitudes  f,  g  and  h. 


...  . .  3f  mrt 
(M+1)3x  +“b~h  =  0 

....  ,3g  mjt , 

(M-'w+Trh=0 

.  .3h  mrt  .. 

MS-2b('+9»'0 


(2.2) 


Defining  the  vector  U  as 


f 

g 

h 


(2.3) 


the  system  (2.2)  can  be  written  in  matrix  form 

Ux  +  A  U  =  0  (2.4a) 

where  the  subscript  x  denotes  partial  derivative  with  respect  to  x  and  where  the  matrix  A 
is 


0 


A  = 


0 


-mrt 
2  Mb 


0 

0 

-mrt 

2Mb 


mrt 

b(M+1) 

mrt 

b(M-i) 

0 


(2.4b) 


Solutions  of  the  form 

U  =  U0e 


(2.5) 


exist  for  p.  equal  to  the  eigenvalues  of  A,  obtained  from  the  solutions  of  det  (A  -  ja.1)  =  0, 
leading  to 


^2,3 


(2.6) 


with  P2  =  1  -  M2. 

The  corresponding  amplitudes  are  the  eigenvectors  of  A,  denoted  by  U('),  and  are 
given  by 


1 

p 

-P 

M+1 

M+1 

-1 

u(2)  = 

-P 

1-M 

u(3)  = 

P 

1-M 

0 

1 

1 

(2.7) 


Hence,  the  general  solution  is  written  as  a  linear  combination  of  the  above  eigenvectors. 


3 

U  =  X  C(U(i)e  M 

i=1 

or  explicitly, 


f 

1 

P 

-p 

M+1 

M+1 

g 

-c3 

-1 

+  C  2 

-P 

1-M 

mn 

e  bpx +c3 

P 

1-M 

h 

0 

1 

1 

(2.8) 


(2.9) 


The  C-coefficients  are  determined  by  matching  these  solutions  with  the  numerical 
solution  at  the  boundary  of  the  computational  domain.  Denoting  by  fo,  go.  ho  the 
numerical  values  of  the  perturbation  wave  amplitudes  at  the  boundary,  taken  as  the  line 
x— 0,  the  C-  coefficients  satisfy  the  following  system: 


C'+MU^-MirC’=(» 

-c'-l£tc’+l&c>-«. 

C2+  C3=  h0 


(2.10) 


The  solutions  are  given  as 


(2.11) 


C»  =  2M,(I+M)fo+(1-M)g0l 

^2=  2^0'  (^o+ 

C3=  ^  |ho+  (f0+  So)^) 

In  order  to  determine  the  boundary  conditions  valid  at  the  finite  distance  location  of  the 
boundaries  of  the  computational  domain,  corresponding  to  uniform  flow  conditions  at 
infinity,  we  have  to  look  at  the  incoming  and  outgoing  characteristics. 

From  the  properties  of  the  characteristic  variables  it  is  known,  see  for  instance  Hirsch 
(1989)  for  a  general  presentation,  that  w'2  and  w’3  are  characteristics  propagating  from 
left  to  right  (for  positive  u),  while  w'4  is  propagating  from  right  to  left  for  a  subsonic 
flow,  since  they  correspond  resp.  to  wave  speeds  u,  u+c  and  u-c. 

Hence,  in  order  to  determine  the  far  field  disturbances  we  have  to  express  that  the 
amplitudes  of  the  incoming  characteristic  perturbations  are  zero  at  infinity,  leading  to  a 
correction  on  the  physical  boundary  conditions  for  finite  distances,  and  that  the 
amplitudes  of  the  outgoing  characteristics  are  defined  by  the  numerical  solution  at  the 
boundary.  These  equations  are  now  examined  separately  for  an  inlet  and  an  exit 
boundary. 


2.1.  Inlet  boundary  conditions 

At  an  inlet  section,  w'2  and  w'3  are  incoming  characteristics,  while  w'4  is  propagating 
from  inside  the  computational  domain  towards  the  boundary,  for  a  subsonic  infiow,  see 
figure  2.2. 

Referring  to  the  solutions  (2.1)  this  implies  that  the  amplitudes  h  and  f  must  decay  to 
zero  for  x-»-oo  and  consequently,  from  equation  (2.9),  that  the  coefficients  Cj  and  C2 
have  to  vanish  since  they  are  associated  resp.  with  a  constant  or  an  exponentially  growing 
amplitude.  Introduced  in  the  relations  (2.1 1),  these  conditions  determine  the  perturbations 
fo  and  ho  in  function  of  go,  which  in  turn  is  obtained  from  the  internal  numerical  solution. 

The  following  relations  are  therefore  valid  at  the  inlet  boundary,  instead  of  fo  =  ho  =  0, 


1-M 

l°  1+M  8° 

.  _  1-M  _  P 

ho-  ~ ~80-  .  ,.go 
P  1+M 


(2.12) 


and  the  amplitude  C3  of  the  remaining  wave  is 


C3=ho-T+Mgo 


(2.13) 


The  first  order  perturbation  solution  in  the  region  upstream  of  the  inlet  boundary  is 
then  completely  defined  for  each  harmonic  as  (x  <  0) 


f 

-P 

M+1 

g 

90P 

P 

~  1+M 

1-M 

h 

1 

nw 

ebpx 


(2.14) 


The  corrected  boundary  treatment  is  therefore  established  as  follows,  assuming  that  the 
incoming  flow  is  isentropic: 

i)  Develop  the  numerically  computed  distribution  of  the  characteristic  perturbation 


-(q-qj 


(2.15) 


in  a  Fourier  series  in  y  along  the  inlet  section.  With  the  definitions  (1.26)  and  the  first 
relation  (2.12),  we  have 


q'  =  u'  = 


1  +  M. 


w4  = 


-1 


1  +  M. 


9o 


(2.16) 


and  hence  it  is  easier  to  develop  the  perturbations  of  the  x-component  of  the  velocity  in 
the  inlet  section,  in  a  cosine-Fourier  series 


u’  =  XAmcos^  (217) 

m  u 

In  practice,  4  to  5  harmonics  are  sufficient  for  the  required  accuracy. 

ii)  At  a  duct  inlet  section,  it  is  customary  to  fix  stagnation  pressure  and  temperature,  as 
well  as  the  inlet  flow  angle.  The  present  treatment  provides  a  correction  to  the  uniform 
inlet  angle  due  to  the  finite  distance  of  the  boundary  and  is  defined  by  applying  the 
Am-coefficients  to  the  following  development,  applying  equation  (2.1 3), 

w2  =  v'-q«®'  =  ~  P  (2.i8) 

m  u 

Note  the  sine-expansion  in  this  formula. 

iii)  The  pressure  disturbances  at  inlet  are  directly  obtained  from  u'  by  application  of 
equation  (1.31)  for  isentropic  conditions,  that  is 

P'  =  *u’P„Uoo  (2.19) 


without  the  necessity  of  performing  a  Fourier  expansion. 


In  summary,  the  numerically  computed  x-component  of  the  velocity 
perturbations  with  respect  to  the  far  field  velocity  IJM  is  expanded  in  a 
cosine-Fourier  series  along  the  inlet  section.  The  resulting  Fourier 
coefficients  are  applied,  after  multiplication  by  (-[}),  as  coefficients 
of  a  sine-Fourier  expansion  in  order  to  obtain  the  flow  angle  disturbances 
along  the  inlet  section. 

The  corrected  flow  angles  are  then  applied  as  new  boundary  conditions. 


2.2.  Exit  boundary  conditions 

At  an  exit  section,  w'2  and  w’3  are  outgoing  characteristics,  while  w'4  is  an  incoming 
characteristic  for  subsonic  exit  velocities,  figure  2.3.  Consequently,  only  one  physical 
boundary  condition  can  be  imposed  at  a  subsonic  exit,  generally  the  static  pressure. 

Referring  to  the  solutions  (2.1)  this  implies  that  the  amplitude  g  must  decay  to  zero  for 

x— >00  and  consequently,  from  equation  (2.9),  that  the  coefficient  C3  has  to  vanish  since  it 
is  associated  with  an  exponentially  growing  amplitude.  Introduced  in  the  relations  (2.1 1), 
this  condition  determines  the  perturbation  go  in  function  of  fyand  ho,  which  in  turn  are 
obtained  from  the  internal  numerical  solution. 

The  following  relations  are  therefore  valid  at  the  exit  boundary,  instead  of  go  =  0, 


9n  =  -rMh0-fc 


(2.20) 


and  the  amplitudes  of  the  remaining  waves  are 


r  f  1"Mh  f  p  h 

C,  ”  '0  h°=  °-  iTM  ho 

C  5  =  h  n 


(2.21) 


I  he  first  order  perturbation  solution  in  the  region  downstream  of  the  exit  boundary  is  then 
completely  defined  for  each  harmonic  as  (x  >  0) 


f 

-3 

1 

Mr-1 

g 

-  ^0  “  h0) 

-1 

ho 

P 

i-M 

h 

0 

1 

nm 

e  8  p  * 


(2.22) 


Note  that  the  coefficient  Cj  is  the  amplitude  of  a  constant  wave  in  x  and  represents  a  non¬ 
decaying  contribution,  in  other  words  a  purely  transported  quantity.  It  is  shown  in  the 
following  that  the  vorticity  disturbances  are  proportional  to  C|. 

The  corrected  boundary  treatment  is  therefore  established  as  follows, 

i)  Develop  the  numerical  distribution  of  the  characteristic  perturbations 


w0  =  v'  =  q  6' 

C.  oo 


(2.23) 


in  a  Fourier  series  in  y  along  the  exit  section.  With  the  definitions  (1.26)  and  the 
relations  (2.21),  the  coefficients  (p™  cM  (g+f)/2]  are  equal  to  the  Fourier  coefficients  B 
of  the  expansion  of  the  pressure  disturbances  in  a  cosine  series  and  are  connected  to 
the  v'-expansion  by  application  of  equation  (2.20).  With 


P’ 


=  lBm 


cos 


mtty 


(2.24a) 


v'  =  q~9  =  X  Amsin^~ 

m  u 

we  have, 

“  Poo^ooy 


(2.24b) 


(2.25) 


Hence,  in  order  to  obtain  the  pressure  corrections  due  to  the  finite  distance  of  the  exit 
boundary,  it  is  sufficient  to  develop  only  the  v'-perturbations  in  a  sine-Fourier  series 


ii)  The  corrections  to  the  uniform  pressure  boundary  condition  are  defined  by  applying 
the  Am-coefficients  to  the  following  development,  leading  to  the  new  boundary 
condition 


P  P 

r  oo~  oo  ^  . 

P  =  P - —  L  A 

00  p  . 


rrmy 


(2.26) 


Note  the  cosine-expansion  in  this  formula. 

In  summary,  the  numerically  computed  y-component  of  the  velocity 
perturbations  is  expanded  in  a  sine-Fourier  series  along  the  exit  section. 
The  resulting  Fourier  coefficients  are  applied,  after  multiplication  by  the 
appropriate  coefficients,  to  a  cosine-Fourier  expansion  in  order  to  obtain 
the  new  pressure  boundary  condition. 

The  solution  (2.22)  shows  that  a  constant  amplitude  wave,  first  term  in  (2.22),  is 
maintained  at  downstream  infinity.  This  wave  is  connected  to  the  C|  coefficient  and  is  due 
to  the  vorticity  waves  which  generate  the  eigenvalue  pi=()  in  equation  (2.6).  This  can  also 
be  shown  from  a  direct  computation  of  the  vorticity  perturbation 
,  3u'  civ' 

0)  =  - v  - 

r)y  dx 


(2.27) 


which  can  be  expressed  in  function  of  the  characteristic  variables  and  the  expansions 
(2.1),  for  an  arbitrary  Fourier  mode  m,  by 


,  ,rrm  f  -  g  5h  .  mny 


(2.28) 


When  applying  the  above  relations  (2.22),  valid  at  the  exit  station,  we  obtain  after  some 
algebraic  manipulations. 
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(2.29) 


This  vorticity  is  generated  in  the  computational  domain  and  is  either  due  to  a  physical 
source,  such  as  a  non-uniform  shock  intensity,  or  to  numerical  dissipation  acting  as 
numerical  vorticity  sources  if  the  flow  is  isentiopic. 

In  both  cases,  the  generated  vorticity  is  transported  out  of  the  domain,  without 
decaying  with  downstream  distance. 

The  value  of  the  Ci  coefficient,  calculated  from  the  above  boundary  procedure,  is  a 
measure  of  the  numerical  accuracy  for  isentropic  flows,  since  it  will  be  a  measure  of  the 
numerically  generated  vorticity. 


2.3.  Applications:  Subsonic  channel  flow 

The  treatment  of  the  far  field  boundary  conditions,  described  above,  is  applied  to  a 
channel  flow  with  subsonic  inlet  and  outlet  conditions.  The  channel  is  formed  by  a  lower 
wall  with  a  sinusoidal  shape,  in  order  to  have  a  smooth  variation  of  slope,  and  a 
rectilinear  upper  wall. 

Calculations  are  performed  in  the  reference  domain,  shown  in  figure  2.4,  with  61 
points  in  the  strearnw  i.se  direction  and  compared  to  results  for  a  restricted  domain  of  31 
streamwise  points,  where  the  boundaries  have  been  taken  much  closer  to  the  central  part 
of  the  channel. 

The  improved  boundary  treatment  is  applied  to  an  Filler  code  developed  at  the  VIJB. 
This  code,  see  l.acor  and  Hirsch  (1988a),  (1988b)  for  a  more  detailed  description,  solves 
the  time  dependent  Euler  equations  in  conservation  form  on  a  stmetured  mesh,  applying  a 
finite  volume,  cell-centered  discretization.  The  fluxes  are  calculated  with  an  upwind 
method  based  on  flux  splitting,  extended  to  second  order  by  variable  extrapolation 
(MUSCL  approach)  and  the  TVD  property  is  maintained  via  the  introduction  of  limiters. 
The  time  integration  is  implicit  and  the  corresponding  implicit  operators  are  solved  by 
relaxation  methods,  coupled  to  a  multigrid  technique. 

Figure  2.5  shows  the  isoMach  number  distribution  in  the  central  part  of  the  channel, 
comparing  the  results  obtained  for  the  extended  and  restricted  domains,  the  latter  with  the 
corrected  boundary  treatment.  Figure  2,5a  is  obtained  with  the  standard  boundary 
conditions,  imposing  uniform  pressure  at  exit  and  zero  flow  angle  at  inlet,  but  only  the 
central  part  is  shown.  When  the  corrected  conditions  are  introduced,  figure  2.5b,  the 
computations  on  the  restricted  domain  become  very  close  to  the  solution  on  the  extended 
domain. 

Figure  2.6  shows  the  Mach  number  distribution  on  the  lower  waii,  comparing  the 
results  obtained  for  the  extended  and  restricted  domains  without  and  with  the  corrected 
boundary  treatment.  Figure  2.6a  is  obtained  with  the  standard  boundary  conditions, 
imposing  uniform  pressure  at  exit  and  zero  flow  angle  at  inlet.  The  errors  introduced  by 
these  conditions  when  the  limits  of  the  computational  domain  arc  too  close  is  clearly  seen. 


When  the  corrected  conditions  are  introduced,  figure  2.6b,  the  computations  on  the 
restricted  domain  become  very  close  to  the  reference  solution  on  the  extended  domain. 

2.3.2.  Transonic  channel  flow 

The  boundary  corrections  defined  in  section  2.2  do  not  require  the  flow  to  be  ientropic 
at  outlet,  and  therefore  the  method  should  apply  also  to  non-isentropic  conditions,  as 
occurs  when  a  shock  appears  in  the  channel.  If  the  inlet  Mach  number  is  increased 
compared  to  the  previous  case,  a  curved  shock  appears,  as  seen  on  figure  2.7,  resulting  in 
a  non  uniform  entropy  downstream  of  the  shock.  Figure  2.7  shows  the  isoMach  number 
distributions  in  the  central  part  of  the  long  channel,  comparing  the  results  obtained  for  the 
extended  and  restricted  domains,  the  latter  with  uncorrected  (b)  and  corrected  (c) 
boundary  treatment. 

The  Mach  number  distributions  on  the  lower  and  upper  walls  are  shown  on  figure  2.8, 
for  the  three  cases  of  figure  2.7.  There  is  a  shift  in  the  shock  position  by  one  mesh  cell, 
which  is  not  very  significant  even  on  this  relatively  coarse  mesh.  The  improvement  due  to 
the  boundary  corrections  is  clearly  seen.  Another  measure  of  the  corrections  concerns  the 
inlet  angles;  the  corrected  inlet  angle  for  the  short  channel  is  2.6  degrees,  to  be  compared 
with  the  value  of  2.7  deg.  calculated  along  the  same  section  of  the  long  channel,  while  in 
the  uncorrected  case  the  inlet  angle  is  fixed  at  zero  degrees.  Another  display  of  the  effects 
of  the  boundary  treatment  is  shown  on  figure  2.9  where  the  Mach  number  profiles  are 
compared  at  inlet  and  outlet  of  the  short  channel.  The  differences  between  the  dashed  lines 
and  the  plus-symbols  indicate  the  amplitude  of  the  corrections  on  the  short  channel,  while 
the  solid  line  is  the  reference  value  from  the  long  channel.  The  small  difference  between 
the  latter  and  the  corrected  values  of  the  short  channel  computation  (+  symbols)  is 
probably  due  to  the  fact  that  the  boundaries  of  the  long  channel  have  not  been  taken  far 
enough. 


3.  FAR  FIELD  BOUNDARY  CONDITIONS  FOR  CASCADE  FLOWS 

The  method  described  in  section  2  is  extended  to  cascade  flows  with  a  geometrical 
configuration  shown  in  figure  3.1.  Cascade  flows  are  characterized  by  the  fact  that  the 
passage  between  two  blades  is  repeated  indefinitely  in  the  direction  of  the  cascade  front, 
that  is  the  y-direction,  with  a  periodicity  equal  to  the  pitch  s  (the  y-direction  corresponds 
to  the  tangential,  wheel  speed,  direction  in  the  turbomachine).  Due  to  this  periodicity,  the 
computational  domain  is  limited  by  the  blade  surfaces  and  two  periodic  extensions.  At 
inlet  the  boundaries  AE  and  BF  are  periodic,  that  is  all  physical  flow  properties  at 
corresponding  points  P  and  P'  are  identical.  Similar  properties  exist  at  exit  along  the 
segments  GC  and  HD.  The  inlet  and  exit  sections  AB  and  CD  have  to  be  parallel  to  the 
cascade  front  in  order  to  satisfy  the  periodicity  conditions. 

The  standard  physical  boundary  conditions  which  are  valid  strictly  at  infinity,  are 
uniform  values  of  flow  angle,  stagnation  pressure  and  temperature  at  inlet  and  static- 
pressure  at  exit. 

The  objective  of  the  present  approach  is  to  deHne  corrections  on  these  conditions  due 
to  the  finite  distance  of  the  limits  of  the  computational  domain.  Since  stationary  flows  are 
being  considered,  stagnation  temperature  remains  constant  and  hence  only  the  flow  angle 
has  to  be  corrected  at  inlet  and  the  static  pressure  at  outlet. 

The  basic  equations  are  given  by  the  linearized  system  ( 1 .24),  written  for  an  arbitrary 
wave  propagation  direction  tc  and  we  look  for  stationary,  perturbation  solutions  of  the 
system  (1.24)  in  the  region  between  the  inlet  (or  exit)  station  of  the 
computational  domain  and  the  boundary  at  infinity. 


Since  entropy  is  purely  convected  and  decoupled  from  the  other  equations,  we  can 
solve  separately  for  the  entropy  perturbation  and  remove  the  corresponding  equation  from 
the  system  (1.24). 

Disturbance  solutions  by  separation  of  variables  are  sought  for  the  remaining  variables 
in  the  external  region,  with  a  complex  Fourier  series  in  y,  such  as  to  satisfy  periodicity 
with  a  period  equal  to  s. 
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Compared  to  the  duct  case,  the  amplitudes  H,  F  and  G  are  complex  quantities,  written  as 
H  =  ht - i h2 

F  =  f1-if2  (3.2) 

G  =  g,-ig2 

where  the  subscript  m  has  been  removed  for  clarity  of  the  notations. 

The  complex  form  of  the  solutions  (3.1)  is  a  compact  way  of  expressing  periodicity, 
but  obviously  only  the  real  part  has  physical  significance.  The  real  part  of  the  solution  is 
given,  for  an  arbitrary  Fourier  mode  m,  by 


(3.1a) 

(3.1b) 

(3.1c) 


w2  =  h ,  cos<}>y  +  h2  sinijjy 

w3  =  f !  coscjty  +  f2  sin0y  (3.3) 

w4  =  g1  cos<{>y  +  g  2  sin<t>y 


where  the  phase  angle  <j>  is  defined  by 


4>  = 


2jtmy 

~s 


(3.4) 


Introducing  these  solutions  in  the  stationary  form  of  equations  (1.24)  we  obtain  the 
following  system,  for  all  Fourier  modes,  as  a  consequence  of  the  linearity  of  the 
perturbation  equations,  introducing  the  axial  (x)  and  tangential  (y)  Mach  numbers, 
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(3.5) 


Ma=  M^cosQ^ 
Mu=  M^sine^ 


and  defining  a  =  sina  and  b  =  cosa. 


Ma|^-|  ~(F+G)  +  i 0MU  H  +  i | 0  (F+G)  =  0 

(Ma+  b)  ^  -  a  ^-  +  i  0  hH  +  i  0  (Mu  +  3)F  =  0  (3.6) 

r)H 

(Mg-  b)  -  a  -*  i  0  bH  +  i  0  (Mu  -  a)F  =  0 
Defining  the  vector  U  as 


U  = 


H 

F+G 

F-G 


(3.7) 


the  system  (3.6)  can  be  written  in  matrix  form,  after  some  rearrangement, 

A  Ux  +  B  U  =  0  (3.8) 

where  the  subscript  x  denotes  partial  derivative  with  respect  to  x  and  the  matrix  A  is 
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-  2a  Ma  b 
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It  is  of  interest  to  note  at  this  point  that  the  (F-tG)  terms  are  the  Fourier  coefficients  of 
2p'/(pooC  oo)  and  that  the  difference  (F-G)  represents  the  expansion  of  the  velocity 
perturbation  along  the  propagation  direction  ie,  that  is  v'-ic.  Note  also  the  appearance  of 
the  axial  and  tangential  Mach  numbers,  which  are  typical  for  cascade  flows. 

Solutions  of  the  form 

U  =  U0e-M*  (3.10) 

are  sought,  where  (.t  is  equal  to  the  eigenvalues  of  (B-Ap),  leading  to 

Pl=i0T  with  t  =  tanOoo  (3.11) 

while  the  two  other  solutions  are  obtained  from  the  roots  of  the  quadratic  equation,  setting 
p=03., 


(3.12) 


X2(l-M2a)+2iXM^lu-(l-M2u)  =  0 
The  solutions  are 


-ixMaMu±P  x 
^23“  2  ^ 

1-Ma 


(3.13) 


with  p2=l-  M2,  where  M  is  the  Mach  number  of  the  far  field.  The  appearance  of  this  term 
indicates  that  the  considered  solutions  become  oscillatory  in  x  at  supersonic  Mach 
numbers,  since  p  becomes  purely  imaginary  in  this  case. 

It  is  of  importance  to  observe  that  these  eigenvalues  are  completely 
independent  of  the  direction  a  of  the  propagation  vector  ic  .  They  depend 
only  on  the  flow  parameters  in  the  far  field,  namely  Mach  number  and  flow  angle  at 
infinity.  This  is  to  be  expected  since  the  linearized  solutions  of  the  Euler  perturbation 
equations  should  have  a  space  dependency  based  on  physical  quantities.  Only  the 
amplitudes  of  these  waves  will  be  function  of  the  selected  propagation  direction. 

The  denominator  in  equation  (3.13)  is  positive  for  subsonic  axial  velocities,  which  is 
always  the  case  in  practice. 

The  amplitudes  of  the  solution  (3.10)  are  the  eigenvectors  of  the  matrix  (B-A|i), 
associated  to  the  three  eigenvalues  (3.1 1),  (3.13).  They  are  easily  obtained  as 


1  Jzl 

l  i+4X 

2  1+^x 
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where  X  is  either  of  the  the  eigenvalues  defined  by  (3.13)  and  E  =  tana. 

Note  that  these  eigenvectors  are  independent  of  the  phase  angle  4>,  but  depend  on  the 
propagation  direction  a.  Two  of  these  directions  are  of  interest,  the  first  one  corresponds 
to  propagation  directions  along  the  free  stream  at  infinity,  that  is  a=0oo  or  E,=x  and  the 

second  one  to  a  propagation  along  the  axial  direction,  that  is  a=0.  The  consequences  of 
these  choices  are  connected  to  the  values  of  the  perturbations  to  the  acoustic  waves,  as 
described  by  the  third  and  fourth  characteristics  (1.20)  and  (1.21).  In  function  of  the 
cartesian  components  u',  v'of  the  velocity  disturbances,  we  have 


w . 
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P~c~> 


(3.15a) 
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-  u'  cosa-  v'  sina 


(3.15b) 


and  in  the  first  case  the  velocity  disturbances  are  equal  to  q',  the  disturbance  of  the 
velocity  amplitude,  while  in  the  second  case,  the  velocity  perturbations  are  restricted  to  the 
axial  component. 

We  will  consider  in  the  following  the  first  choice,  that  is  a=0oo  or  ^=T,  observing  that 


1+iTp 
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(3.16) 


the  eigenvectors  (3.14)  become, 
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it 

-  M 
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(3.17) 


and  the  general  solution  is  a  linear  combination  of  these  eigenvectors 
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with 


K  = 
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1  -m\ 


L  = 


1  -  M3 


(3.19) 


The  complex  C-cocfficients  are  determined  by  matching  these  solutions  with  the 
numerical  solution  at  the  boundary  of  the  computational  domain.  Denoting  by  F(),  Go, 
Ho  the  numerical  values  of  the  perturbation  wave  amplitudes  at  the  boundary,  taken  as  the 
line  x=0,  we  define  real  and  imaginary  parts  of  the  amplitudes 


^0  =  ^01  -*^02 

Po  =  F0+ G0  =  p0, -ip02  (3.20) 

Q0  =  Po-^0  =  Cloi~'^02 
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The  coefficients  pot  and  po2  are  the  Fourier  coefficients  of  the  pressure  fluctuations, 
while  qoi  and  qo2  correspond  to  the  expansion  of  the  velocity  disturbance  q',  as  can  be 
seen  from  equations  (3. 15)  above. 

The  C-  coefficients  satisfy  the  following  system, 


P(C3-C2)  =  H0 
-  M  (C2+  C3 )  =  P0 

c1  +  c2  +  c3=q0 


(3.21) 


Identifying  real  and  imaginary  parts  and  solving  for  the  unknowns  leads  to  the  relations 
for  the  real  values,  marked  by  a  superscript  (R). 
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(3.22a) 


MCr=Mq01  +  p01 

and  for  the  imaginary  components,  marked  by  a  superscript  (I), 
hoi  .  P02 


2Cf=-2 
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(3.22b) 


Referring  to  equation  (3.3),  it  is  seen  that  the  real  parts  of  Po  and  Qq  correspond  to  the 
coefficients  of  the  cosine  teams  in  the  Fourier  expansion,  while  I102  is  the  coefficient  of 
the  sine  terms  in  the  expansion  of  w'2.  That  is 


w2  =  =  Z  (hoi  cos<}>y  +  h02  sin<t>y) 

m 

Q  =  2q'  =  X  (Qoi cos^  +  Q02  sin<J>y) 


(3.23) 


In  order  to  determine  the  boundary  conditions  valid  at  the  finite  distance  location  of  the 
boundaries  of  the  computational  domain,  corresponding  to  unifomi  flow  conditions  at 
infinity,  we  proceed  in  the  same  way  as  in  section  2  and  express  that  the  amplitudes  of  the 
incoming  characteristic  perturbations  are  7cro  at  infinity,  leading  to  a  correction  on  the 
physical  boundary  conditions  for  finite  distances,  and  that  the  amplitudes  of  the  outgoing 


characteristics  are  defined  by  the  numerical  solution  at  the  boundary.  These  equations  are 
now  examined  separately  for  an  inlet  and  an  exit  boundary. 


3.1.  Inlet  boundary  conditions 

At  an  inlet  section,  w‘2  and  w'3  are  incoming  characteristics,  while  w'4  is  propagating 
from  inside  the  computational  domain  towards  the  boundary,  for  a  subsonic  inflow,  see 
figure  3  2. 

Referring  to  the  solutions  (3.1)  this  implies  that  the  amplitudes  H  and  F  must  decay  to 
zero  for  x-*-«>  and  consequently  that  the  coefficients  Cj  and  C2  have  to  vanish  since 
they  are  associated  resp.  with  a  constant  or  an  exponentially  growing  amplitude. 
Introduced  in  the  relations  (3.22),  these  conditions  determine  the  perturbations  Fo  and  Ho 
in  function  of  Go,  which  in  tum  is  obtained  from  the  internal  numerical  solution. 

The  following  relations  are  therefore  valid  at  the  inlet  boundary,  instead  of  Fo  =  Ho  = 

0, 
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(3.24) 


The  first  order  perturbation  solution  in  the  region  upstream  of  the  inlet  boundary  is  then 
completely  defined  for  each  harmonic  as  (x  <  0) 


u  =  (q0i-icin?) 
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(3.25) 
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The  perturbations  of  the  outgoing  characteristic  w'4  can  be  expressed  in  function  of  the 
disturbances  of  the  magnitude  of  the  velocity  q'  by  applying  the  above  relations  (3.24), 
leading  to 

w4=|  (P  — Q)  =  — 1(1  +  MJCU-  (1  +  MJq'  (3.26) 

Therefore,  expanding  the  fourth  characteristic  is  equivalent  to  a  Fourier  expansion  of  the 
velocity  disturbances  q'. 

The  corrected  boundary  treatment  is  therefore  established  as  follows,  assuming  that  the 
incoming  flow  is  isentropic: 

i)  Develop  the  numerically  computed  distribution  of  the  perturbations  of  the  velocity 
magnitude  in  the  inlet  section,  in  a  complete  Fourier  series 

q'  =  X(Amcos<}>y  +  Bmsin<!>y)  (3.26) 


In  practice,  4  to  5  harmonics  are  sufficient  for  the  required  accuracy. 


ii)  At  a  cascade  inlet  section,  it  is  customary  to  fix  stagnation  pressure  and  temperature, 
as  well  as  the  inlet  flow  angle.  The  present  treatment  provides  a  correction  to  the 
uniform  inlet  angle  due  to  the  finite  distance  of  the  boundary  and  is  defined  by 
applying  the  Am  and  Bm-coefficients  to  the  following  development,  applying  equation 
(3.23), 

q^'  =  pX(BmCOS4>y-Amsin<{>y)  (3-27) 

m 

iii)  The  pressure  disturbances  at  inlet  are  directly  obtained  from  q'  by  application  of 
equation  (1.31)  for  isentropic  conditions,  that  is 

P'  =  -q,PootU  (3-28) 

without  the  necessity  of  performing  a  Fourier  expansion. 

In  summary,  the  numerically  computed  perturbations  of  the  velocity 
magnitude  are  expanded  in  a  complete  Fourier  series  along  the  inlet 
section.  The  resulting  Fourier  coefficients  are  applied,  after  multiplication 

by  p,  to  a  modified  Fourier  expansion  in  (3.27)  in  order  to  obtain  the 
flow  angle  disturbances  along  the  inlet  section. 

The  corrected  flow  angles  are  then  applied  as  new  boundary  conditions. 


2.2.  Exit  Imuntiarv  conditions 

At  an  exit  section,  w'2  and  w'3  are  outgoing  characteristics,  while  w'4  is  an  incoming 
characteristic  for  subsonic  exit  velocities.  Consequently,  only  one  physical  boundary 
condition  can  be  imposed  at  a  subsonic  exit,  generally  the  static  pressure. 

Referring  to  the  solutions  (3.1)  this  implies  that  the  amplitude  G  must  decay  to  zero  for 

x— and  consequently,  from  equation  (3.1  8),  that  the  coefficient  C3  has  to  vanish  since 
it  is  associated  with  an  exponentially  growing  amplitude.  Introduced  in  the  relations 
(3.21),  this  condition  determines  the  perturbation  Go  in  function  of  Fo  and  Ho,  which  in 
turn  are  obtained  from  the  internal  numerical  solution. 

The  following  relations  are  therefore  valid  at  the  exit  boundary,  for  the  amplitudes  of  the 
remaining  waves,  instead  of  Go  =  0, 


Ci  =  Oo+My 


C2  =  - 


2|Ho 

~P 


(3.29) 


The  first  order  perturbation  solution  in  the  region  downstream  of  the  exit  boundary  is  then 
completely  defined  for  each  harmonic  as  (x  >  0) 


U  =  (Qo  +  w 
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(3.30) 


Note  that  the  coefficient  Cj  is  the  amplitude  of  a  sinusoidal  wave  in  x  and  represents  a 
non-dccaying  contribution,  in  other  words  a  purely  transported  quantity.  As  shown  in 
section  2,  the  vorticity  disturbances  are  proportional  to  Cj. 

The  corrected  boundary  treatment  is  therefore  established  as  follows, 


i)  Develop  the  numerical  distribution  of  the  flow  angle  perturbations  in  a  Fourier  series 
in  y  along  the  exit  section. 


q~0'  =  X  (AmC0S<j>y  +  Bmsin<j>y) 


(3.31) 


With  the  definitions  (3.23)  and  the  relations  (3.29)  the  expansion  of  the  pressure 
disturbances  is  directly  obtained  from 


p  =  2^— -  =  X  (p0i  cos<i>y  +  p02  sin<t>y) 


(3.32) 


ii)  The  corrections  to  the  uniform  pressure  boundary  condition  are  defined  by  applying 
the  A^,  Bm  coefficients  to  the  following  development,  leading  to  the  new  boundary 
condition 


P  P 

1  oo  1  c. 


p  =  p  - 


p 


X(Bmcos4>y-  Amsin<j>y) 


(3.33) 


In  summary,  the  numerically  computed  flow  angle  perturbations  are 
expanded  in  a  complete  Fourier  series  along  the  exit  section.  The 
resulting  Fourier  coefficients  are  applied,  after  multiplication  by  the 
appropriate  coefficients,  to  a  Fourier  expansion  in  order  to  obtain  the  new 
pressure  boundary  condition. 


3.3.  Ann  I  i  cations  to  cascade  flows 


The  treatment  of  the  far  field  boundary  conditions,  described  above,  is  applied  to  a 
cascade  with  subsonic  inlet  and  outlet  conditions.  The  blade  is  formed  by  a  suction 
surface  with  a  sinusoidal  shape,  in  order  to  have  a  smooth  variation  of  slope,  and  a 
rectilinear  pressure  surface.  The  shape  is  identical  to  the  one  used  in  the  duct 
computations  of  section  2. 

Calculations  are  performed  in  the  reference  domain,  shown  in  figure  3.3a  with  91 
points  in  the  streamwise  direction  and  compared  to  results  for  a  restricted  domain  of  55 
streamwise  points,  figure  3.3b,  where  the  boundaries  have  been  taken  much  closer  to  the 
central  part  of  the  channel.  The  blade  extends  from  -0.5  to  0.5  and  the  long  channel 
extends  from  -1  to  +1;  the  incidence  flow  angle  is  taken  equal  to  5  degrees. 


Figure  3.4  shows  the  isoMach  number  distribution  in  the  central  part  of  the  channel, 
comparing  the  results  obtained  for  the  extended  and  restricted  domains,  the  latter  with  the 
corrected  boundary  treatment.  Figure  3.4a  is  obtained  with  the  standard  boundary 
conditions  on  the  extended  computational  region,  imposing  uniform  pressure  at  exit  and  5 
degrees  flow  angle  at  inlet,  but  only  the  central  part  is  shown.  Figure  3.4b  is  obtained 
with  the  standard  boundary  conditions  on  the  restricted  computational  region,  imposing 
the  same  uniform  boundary  conditions.  When  the  corrected  conditions  are  introduced, 
figure  3.4c,  the  computations  on  the  restricted  domain  become  very  close  to  the  solution 
on  the  extended  domain.  This  can  be  noticed  in  particular  by  the  improved  symmetry  of 
the  distribution  in  the  second  case.  Another  effect  appears  on  the  exit  flow  angles  which  is 
-5.5  for  the  non  corrected  short  domain  and  -5. 1  after  corrections. 

Figure  3.5  shows  an  enlarged  view  of  the  Mach  number  distribution  on  the  lower 
(suction  side)  boundary  comparing  the  results  obtained  for  the  extended  and  restricted 
domains  without  and  with  the  corrected  boundary  treatment.  The  views  cover  the  region 
between  the  leading  and  trailing  edges  resp.  and  the  limits  of  the  computational  domain. 
Figure  3.5a  shows  the  inlet  region,  while  figure  3.5b  shows  the  exit  region.  The  solid 
line  is  the  reference  computation  with  the  extended  region,  the  dotted  line  refers  to  the 
short  channel  and  uniform  boundary  conditions,  while  the  crosses  are  obtained  from  the 
short  channel  with  the  corrected  boundary  treatment.  The  improvement  is  clearly  seen. 

An  other  illustration  of  the  validity  of  the  theory  can  be  seen  from  Figure  3.6,  where 
the  Mach  number  profiles  along  the  inlet  and  exit  station  of  the  short  domains  are  shown 
for  the  same  three  cases.  The  improvement  is  indeed  spectacular. 


Conclusions 

A  method,  based  on  the  linearized  solutions  of  the  full  system  of  Euler  equations  has  been 
developed,  for  channel  and  cascade  flows,  in  order  to  correct  the  uniform  boundary 
conditions,  strictly  valid  at  infinity,  for  the  finite  distance  of  the  limits  of  the 
computational  domain.  The  linearized  solutions  are  obtained  by  separation  of  variables 
with  a  Fourier  expansion  in  the  coordinate  along  the  inlet  and  exit  stations  and  an 
exponential  variation  in  the  axial  direction. 

The  corrections  of  the  boundary  conditions  are  derived  on  the  basis  of  characteristic 
theory,  expressing  that  the  incoming  characteristic  disturbances  have  to  vanish  at  infinity. 
The  outgoing  characteristic  disturbances  are  obtained  from  a  Fourier  expansion  of  the 
numerical  solution  at  the  boundaries.  The  Fourier  coefficients  obtained  from  this 
expansion  are  used  to  generate  the  non  uniform  corrections  on  the  physical  boundary 
conditions  such  as  flow  angle  or  pressure. 

Computations  on  duct  and  cascade  flows  show  the  correctness  and  accuracy  of  the 
method,  for  isentropic  and  non-isentropic  conditions,  allowing  a  considerable  reduction 
of  the  size  of  the  computational  domain. 

Some  further  work  should  be  done  to  extend  the  approach  to  cases  where  shocks  cross 
the  boundaries,  in  order  to  cover  most  of  the  cases  occurring  in  practice. 

acknowledgment 

This  work  was  initiated  and  partially  supported  by  Naval  Air  Systems  Command  (G. 
Derderian-  Program  Manager)  under  contracts  N62271-86-M-0202  and  N6227I  87-M- 
0215,  from  the  Naval  Postgraduate  School,  Monterey,  Ca,  where  the  technical  monitor 
was  Prof.  R.P.Shreeve. 

The  author  would  also  like  to  thank  Mr.  F.  Alcrudo,  Ph  I).  student  at  the  VUB,  for  his 
help  in  performing  the  computations  reported  in  sections  2  and  3. 


References 


BAYLISS  A.,  TURKEL  E.  (1982).  "Far  Field  Boundary  Conditions  for  Compressible 
Flows".  Journal  Computational  Physics,  Vol.48,  pp.  182- 199. 

HIRSCH  Ch.,  LACOR  C.,  DECONINCK  H.  (1987).  "Convection  Algorithms  Based  on 
a  Diagonalization  Procedure  for  the  Multidimensional  Euler  Equations"  Proc.  AIAA  8th 
Computational  Fluid  Dynamics  Conference,  AIAA  Paper  87-1163,  pp.667-676. 

HIRSCH  Ch.  (1988).  Numerical  Computation  of  Internal  and  External  Flows.  Vol.l: 
Fundamentals  of  Numerical  Discretization.  J.  Wiley  ,  London 

HIRSCH  Ch.  (1989).  Numerical  Computation  of  Internal  and  External  Flows.  Vol.2: 
Computational  Methods  for  Inviscid  and  Viscous  Flows.  J.  Wiley  ,  London 

LACOR  C.,  HIRSCH  Ch.  (1988a).  "  3D  Computations  of  Complex  Flow  Systems". 
Paper  presented  at  the  16th  1CAS  Conference,  Jerusalem. 

LACOR  C.,  HIRSCH  Ch.  (1988b).  "Numerical  Simulation  of  the  Three-Dimensional 
Flow  around  a  Butterfly  Valve”.  In  Fluid  Dynamics  in  Non-Rotating  Turbomachinerv 
Components.  ASME  Winter  Annual  Meeting,  Proc.  FED,  Vol.  69,  pp.  1- 1 2. 

ROE  P.L.  (1986).  "Remote  Boundary  Conditions  for  Unsteady  Multidimensional 
Aerodynamic  Computations".  NASA  CR-17821 1, 1CASE  Report  86-75,  NASA  Langley 
Research  Center,  Hampton, Va. 

THOMAS  J.L.,  SALAS  M.D.  (1986).  "Far  Field  Boundary  Conditions  for  Transonic 
Lifting  Solutions  to  the  Euler  Equations".  AIAA  Journal  Vol. 24,  pp.  1074- 1080. 

VERHOFF,A.(1985).  "Modeling  of  Computational  and  Solid  Surface  Boundary 
Conditions  for  Fluid  Dynamics  Calculations".  AIAA  Paper  85-1496  CP,  AIAA  7th 
Computational  Fluid  Dynamics  Conference. 

VERHOFF,A.(1988)."Far  Field  Computational  Boundary  Conditions  for  Internal  Flow 
Problems".  Naval  Postgraduate  School  Report,  NPS  67-88-001CR,  Sept  1988. 

VERHOFF  A.,  O'NEIL  P. J.(  1 984).  "A  Natural  Formulation  for  Numerical  Solutions  of 
the  Euler  Equations".  AIAA  Paper  84-0163. 


Wave  propagation  direction 


— * 


Figure  1.1:  Wave  front  surfaces  and  wave  propagation  direction 


Fig.  1.2:  Velocity  variations  in  streamline  coordinates 
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Figure  2. 1 :  Schematic  representation  of  duct  geometry 
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Fig.  2.3:Fxit  region  and  associated  characteristics 
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MESH  DIMENSION  :  SI*  21 

Figure  2.4,  :  Long  (61x21)  and  short  (31x21)  mesh  for  sinusoidal  channel 


Figure  2.7  :  Comparison  of  iso-Mach  lines  in  the  central  part  of  the  sinusoidal 
channel  at  transonic  conditions  :  (a)  long  channel 

(b)  short  channel  and  uncorrected  boundary  conditions 

(c)  short  channel  with  corrected  boundary  conditions 


Figure  2.8  :  Comparison  of  Mach  number  distributions  in  the  central  part  of  the 
sinusoidal  channel  at  transonic  conditions  on  lower  (a)  and  unner 
(b)  walls 

Solid  line  :  long  channel 

Dashed  line  :  short  channel 'and  uncorrected  boundary  conditions 
♦♦symbols  :  short  channel  with  corrected  boundary  conditions 
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Figure  2.9  :  Comparison  of  Mach  number  distributions  along  the  inlet  (a)  and 
exit  (b)  stations  of  the  short  sinusoidal  channel  at  transonic 
conditions 

Solid  line  :  long  channel 

Dashed  line  :  short  channel  and  uncorrected  boundary  conditions 
++  symbols  :  short  channel  with  corrected  boundary  conditions 
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Fig.  3. 1 :  Cascade  geometrical  configuration 


Fig.  3.2:  Cascade  inlet  flow  conditions 
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Figure  3.3  :  Cascade  geometry  and  mesh  for  the  long  and  short  cascade  passage 


J9 


Suet  ion  tide 


/ 


/ 


/■ 

/ 


(a) 


6urt i on  side 


l.«0  1.60  I.W  2  "0  7  .20  2.41)  2  60  2  60  10  1 

>-coiK<Jin4ie 


¥  igure 
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Enlarged  view  of  the  Mach  number  distributions  on  the  suction  side 
boundary  for  the  three  cases  of  figure  3.4 

(a)  inlet  region  be  tween -x=0 . 2 5  and  -0.15 

(b)  exit  region  between  x-0.15  and  0.25 
Solid  line  :  long  channel 

Dashed  line  :  short  channel  and  uncorrected  boundary  conditions 
++  symbols  :  short  channel  with  corrected  boundary  conditions 
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